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ABSTRACT 

Gravitational coupling between a gaseous disk and an orbiting perturber leads to angular momentum 
exchange between them which can result in gap opening by planets in protoplanetary disks and clearing 
of gas by binary supermassive black holes (SMBHs) embedded in accretion disks. Understanding the 
co-evolution of the disk and the orbit of the perturber in these circumstances requires knowledge 
of the spatial distribution of the torque exerted by the latter on a highly nonuniform disk. Here 
we explore disk-satellite interaction in disks with gaps in linear approximation both in Fourier and 
in physical space, explicitly incorporating the disk non-uniformity in the fluid equations. Density 
gradients strongly displace the positions of Lindblad resonances in the disk (which often occur at 
multiple locations), and the waveforms of modes excited close to the gap edge get modified compared 
to the uniform disk case. The spatial distribution of the excitation torque density is found to be 
quite different from the existing prescriptions: most of the torque is exerted in a rather narrow region 
near the gap edge where Lindblad resonances accumulate, followed by an exponential fall-off with 
the distance from the perturber. Despite these differences, for a given gap profile the full integrated 
torque exerted on the disk agrees with the conventional uniform disk theory prediction at the level of 
~ 10%. The nonlinearity of the density wave excited by the perturber is shown to decrease as the wave 
travels out of the gap, slowing down its nonlinear evolution and damping. Our results suggest that gap 
opening in protoplanetary disks and gas clearing around SMBH binaries can be more efficient than 
the existing theories predict. They pave the way for self-consistent calculations of the gap structure 
and the orbital evolution of the perturber using accurate prescription for the torque density behavior. 

Subject headings: accretion, accretion disks — instabilities — (stars:) planetary systems: protoplan- 
etary disks — (galaxies:) quasars: general 
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1. INTRODUCTION. 

Gravitational coupling between a gaseous disk and an 
external orbiting body plays an important role in many 
astrophysical systems, including disk-planet interaction 
in protoplanetary disks, orbital evolution of supermassive 
black hole (SMBH) binaries surrounded by gaseous disks, 
dynamics of accretion disks in cataclysmic variables, and 
so on. This interaction leads to angular momentum ex- 
change between the disk and the perturber, causing the 
orbit of the latter to evolve — an effect responsi ble for 
the p lanet migration in protoplanetary disks (Ward 1986, 
I1997T ) and the black hole inspiral in the case of SMBH 
binaries embedded in gaseous disks (jlvanov et al.l 119991 : 
IGould fc Ri^l200(l 

Tidal disk-satellit43 coupling also modifies the distri- 
bution of the gas surface densit y and can re sult in gap 
opening in protoplanetary disks (Ward 1997) and cavity 
form ation around the SMBH binaries ( Armitag e et al.l 
2002) if the perturber is massive enough. Gap open- 
ing necessarily weakens the tidal coupling between the 
perturber and the disk, and slows the orbital migration 
of the p erturber, sw itching it to the so-called Type II 
regime, (Ward 1986]). Thus, gap opening can be an im- 
portant "parking mechanism" preventing massive plan- 
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ets from migrating all the way to the central star. 

In steady state the gap structure is determined by the 
local balance between the planetary torque acting on the 
disk and the divergence of the external stress T rtp , which 
can be due to magnetorotational instability or some other 
mechanism of intrinsic angular momentum transport in 
the disk: 



dT 
dr 



dr 



(1) 



Here dT / dr\d is the torque density (the amount of torque 
per unit radial distance) deposited in the disk by the 
density waves originally launched by the planetary tide 
closer to the planet. The deposition of the wave angu- 
lar momentum in the disk occurs b y virtue of some linear 
( Takcuchi et al .11 19961 ) or nonlinear (|Goodman fc Rafikovl 
120011 ) damping process. In general the spatial distribu- 
tion of dT/dr\d can be described as 



dT 
dr 



dr 



(2) 



where the intrinsically non-local operator £ describes the 
wave damping, and dT '/dr is the excitation torque density 
- the rate at which planetary potential adds angular 
momentum to the propagating density waves per unit 
radial distance. 

Thus, in order to understand the details of gap 
opening one must understand both the excitation of 
waves and their damping, a key point first highlighted 
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bv iLunine fc Stevensonl (119821) , iGreenberd ([1983), and 
IGoldreich fc Nicholson! (J1989D . Once the description of 
the wave damping (i.e. the explicit form of the operator 
C) is available, the steady state gap structure is fully de- 
termined by the spatial structure of the excitation torque 
density dT j dr. 

In this work, following the approach of Rafikov & 
Petrovich (2012; hereafter RP12) we concentrate only on 
the pr ocess of wave excitat i on, w hich at least in some 
cases (|Goodman fc Rafikovl 1200 1| ) can be studied sep- 
arately from the wave dissipation (when the latter is 
weak) . Our goal is to provide a self-consistent calculation 
of the excitation torque density dT/dr in non-uniform 
disks affected by gap formation. Prior to this work the 
majority of existing studies of the gap and cavity opening 
in disks have adopted the following simple prescription 
for the excitation torque density behavior in non-uniform 
disks: 



dT(r) _ E(r) dT(r) 



dr 



dr 



(3) 



where dT / dr\ u is the excitation torque density in a uni 
form disk with surface density So, and E(r) is the spa- 
tially varying surface density in a non-uniform disk. 
The behavior of dT/dr\ u has bee n previously de 
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rived from direct numerical simulations (Ba te et alj 
de Val-Borro et alJl2006tlD'Angelo fc Lub ow 2008,12010; 
Dong et aLl l2011al fbf) and from analytical linear stud- 
ies (GT80; Muto & Inutsuka 2009; RP12). In partic- 
ular, in_ttien_pi_ojieermg_^udv_of disk-satellite interac- 
tion IGoldreich fc Trema inc (1983) have shown that far 
from the perturber, at radial separations from its orbit 
\ r ~ r p\ ( r p i s the semimajor axis of the circular orbit of 
the planet) exceeding the disk scale height h, the uniform 
disk excitation torque densitjQ dT / dr\ u is given by 



dT 
dr 



->sign(r -r p )C G T80 



(GM p ) 2 Z 



Q 2 \r 
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GT80 
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81 



2Kn - 



Ml 



2.5, (4) 



where M p is the planetary mass, So is the disk surface 
density assumed to be uniform on scales ~ | r — r p | , fl is 
the angular frequency of the disk at r p , and K n is the 
modified Bessel function of order n. Note that strictly 
speaking dT/dr in equation (4) represents the force den- 
sity and not the angular momentum density as it lacks an 
additional factor of r p . This is because our subsequent 
calculations are cast in the shearing sheet framework in 
which the concept of r p is not well defined. 

The general sc aling dT/dr\„. oc \ r — r j" 4 has been 
also confirmed by lLin fc Papaloizo u (1979) using impulse 
approximat ion. A number of studies of gap opening 
by planets (IL iri & Papaloizou 1986; Tril ling et alJll998t 



Brvden et al. 1999; Armitagc ct al. 2002; Varniere et al. . 
20041: iCrida et al.H2006D and orbital evolut ion of SMBH 



binar i es surrounded by gaseous disks dGould fc RbJ 
120001 lArmitage fc Natarajanl 120021 : lLodato et all 120091 : 

4 We assume that GT80 calculation refers to dT/dr and not to 
dT/dr\d since no explicit dissipation mechanism was mentioned in 
their work. 



lAlexander et al.l 1201 1[ ) have used the prescription^ 
([3]) combined with asymptotic GT80 formula ((4} for 
dT/dr\ u , in some cases with a value of the constant 
pre-factor different from Cot so ()Papaloizou fc Lirilll984t 
lArmitage fc Natar aian 200 2). 

Recently, iDong et al.1 (|2011a| ) carried out high- 
resolution numerical simulations of disk-planet interac- 
tion in the two-dimensional shearing sheet geometry. 
Contrary to the prediction Q of GT80 they found that 
the one-sided excitation torque density dT / dr\ u does not 
maintain a fixed sign but rather changes from positive 
to negative at a distance of ~ 3.2/i from the perturber 
- a result recently confirm ed by global simulations of 
iDuffell fc MacFadvenl (|2012D . This negative torque den- 
sity phenomenon was then explained analytically in the 
framework of linear theory by RP12, who traced its ori- 
gin to the overlap of Lindblad resonances in the vicinity 
of the perturber's orbit (which was ignored in GT80). 

This unexpected finding casts serious doubt on the va- 
lidity of simple prescription ([3]) in non-uniform disks. In- 
deed, directly substituting dT/dr\ u computed by RP12 
in equation ([3]) can result in negative total torque ex- 
erted on the disk by the perturber if the gap is suffi- 
ciently wide and deep (in view of RP12 results this is 
certainly true in the extreme case of the gas density be- 
ing exactly zero at separations less than 3.2h from the 
planetary orbit). This conclusion is unphysical (planet 
would not open a gap in the first place), suggesting that 
the oversimplified prescription ([3]) incorporating the disk 
non-uniformity only through the direct multiplication by 
the local surface density (so that density gradients do not 
affect disk-satellite coupling) needs to be revised. 

In this work, we provide a fully self-consistent linear 
calculation of the disk-satellite interaction by explicitly 
incorporating the disk non-uniformity in the fluid equa- 
tions and properly accounting for the effect of density 
gradients on the waveforms of perturbed fluid variables. 
In this way we are able to compute the torque density and 
angular momentum flux in Fourier and physical space 
and demonstrate significant differences with the results 
obtained using the prescription (j3|) , especially at the gap 
edges where the density gradients are large. 

Our paper is structured as follows. In [21 we describe 
the problem setup, in particular the governing equations, 
the assumed gap density profile, and the Lindblad reso- 
nances in non-uniform disks. Our numerical procedure 
and the results for the density wave behavior are de- 
scribed in SJ3J Calculations of the torque exerted by the 
perturber on the disk in Fourier and real space are de- 
scribed in iJ5]and|ni respectively. Finally, we discuss our 
results, including their astrophysical applications, in Sj8] 

2. PROBLEM SETUP. 

We start by deriving a system of linearized equations 
and describing the adopted underlying surface density 
distribution — the ingredients needed to calculate the 
spatial behavior of the perturbed quantities. 

2.1. Basic equations. 

5 It should also be pointed out that in all these studies the ex- 
citation torque density dT/dr was identified with the deposition 
torque density dT/dr \ ^ , thus ignoring the process of wave dissipa- 
tion altogether, see equation {2}. 
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We study the tidal coupling of a planet with 
a non-uniform disk in the sh earing sheet geometry 
(|Goldreich fc Lvnden-Belllll965f) . which allows us to ne- 
glect geometric curvature effects while preserving the 
main qualitative features of the system. We also neglect 
the vertical dimension and assume the disk to be two- 
dimensional. The dynamics of fluid is then governed by 
the following equation^ (e.g. Narayan et al. 1987; here- 
after NGG): 

<9v VP 

-£- + (v • V) v -I- 2Q. x v + 4Aflxe x = — — - V$ p (5) 



9E 

~dt 



V-(vE) = 0, 



(6) 



where v = (u, v) is the fluid velocity with components 
in the x = r — r p (radial) and y (azimuthal) directions 
correspondingly, P is gas pressure, $ p is the planetary 
potential, fl p = fl p e z is the Keplerian rotation frequency 
at the planet position (r = r p ; we will subsequently drop 
the subscript "p" for brevity), and A = (r/2)(dfl/dr) is 
the shear rate at the same location. Additionally, in our 
notation c s is the sound speed and B = CI + A is Oort's 
B constant. 

For an isothermal disk with no planet present (4> = 
0) and for an arbitrary background density profile E = 
Eo(x), the Equations (O and (|6|) have the following exact 
steady-state solution 



Mo = 0, vq = 2Ax + 



_c|_91nEo 
2~n dx 



(7) 



After having specified our steady-state solution 
{T,q,uo,vq}, we introduce the planetary potential as a 
perturbation and study the linear behavior of the per- 
turbed density Ei and velocity field 5v — (ui,v\). By 
ignoring the quadratic perturbed quantities that appear 
in the Equations (O and (JB]) one can get the following 
general system 

du\ ( . c 2 9 In E \ du\ 

-m + ( 2Ax ' 



— - - 2fi«i 
2ft dx J dy 



„ 2 d(Si/So) _ d% 



dx 



dx 



dt'i ( c 2 s d In E \ dvx 

~m + K 2 + 2n^x~ ) ~dy~ + 

op c g a 2 in sp ^ 2 a(Si/s ) a^ Q . 

2B + 2n^x^J Ul + Cs ^dy— = -^ 9) 
9(Si/Eq) | f 2M | c 2 dlnE \ d^/Sp) 
dt \ 2fl dx J dy 

cUnEo ,du 1 + dv 1=() (1Q) 
dx dy 



u\- 



dx 



We then represent all perturbed fluid vari- 
ables via Fourier integrals as {Ei,Ui,Vi,$ p } = 
(2ir)~ l j exp(ik y y){SJ^,u,v,(j)}, which takes care of 
the y-dependence leaving us with the following set of 

6 Analogous equation (3) in RP12 has a typo: the term 4Aflze x 
should be 4AQxe x . This typo does not affect any other part of the 
paper. 



equations in the ^-coordinate only. Thus, by defining 

(11) 
B B !- — V ' '.".T" (12) 



, . „ , , c 2 k v 91nEo 

a(x) = -2Ak y x - ' ' 



^ 2fl 

c 2 <9 2 lnE 



dx 



AQ dx 2 
we get the linearized system in steady state (d/dt = 0) 

2 <9(<SE/E ) d<t> 



2flv 



2Bu 



-la 



5T, <91nE 



So 



dx 



dx 
du 

— + IKyV ■ 

dx 



dx 

-ik y (j) 



A). 



(13) 
(14) 
(15) 



where the Fourier component of the planetary potential 
produced by a point mass at the origin and its derivative 
are given by 



(f>(ky,x)- 



GAL 



K (\kyx\ 



^sgn^MMIM) 



(16) 
(17) 



Solving (fT5|) - (fT5j) for the azimuthal velocity perturba- 
tion v, we obtain the second order ordinary differential 
equation 



o v 
dx 2 

+v 



<%dln(E /A) 



dx dx 
a 2 ft A _^ k y aAdln(aA/A) 
72" "' 



c 2 B 



2B 



dx 



dcj)2B 
dx c 2 

where we have defined 



kyd fc 2 A91n(A/A) 



2B 



dx 



A = l- 



ky(J 



dB ~ dlnE 
dx dx 



-c 2 k 2 A + 4B 2 . 



(18) 



(19) 
(20) 



The perturbed radial velocity u and the surface density 
perturbation <5E are directly expressed in terms of v via 
the following relations: 



£E 



i 

~A 
Eo 
A 



- r 2^_ 

" v s dx 



,dv 

dx 

~ dv 
2B— + kyffAv 
dx 



2Bav + 2Bk„ 



KH 



(21) 
(22) 



One can easily check that in the limit dT,o/dx —t 
equation (fT5)) reduces to equation (5) of RP12 which was 
derived for a uniform disk. 

2.2. Bisk models 

We model the density structure of an axisymmetric gap 
around the satellite orbit by a particular three parameter 
density profile (symmetric with respect to x = 0) 



Eo(x) 



= 1- 



1 



1 + (x/5) n exp - (2.55 /x) 



(23) 
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Fig. 1. — Density profiles (left panel) and epicyclic fr equ encies 
re (right panel) for the gap models described by Eq. I|23| l with 
parameters <5 = 5h, S m i n = 3 X 10 — 4 , and n = 2, 4, 6. 

where the surface density is normalized to its value at 
infinity E^. In this profile the dimensionless parameter 
E m i n controls the depth of the gap, 5 regulates its width, 
while an even positive integer n sets the steepness of the 
density gradient at the gap edge. The density variation 
between EminEoo inside the gap and Soo outside of it is 
predominantly due to the power law dependence on x in 
the denominator of the second term in the right hand side 
of the equation (|23j). The exponential dependence in the 
denominator is introduced only to ensure a flat bottom at 
x < (5 and allows us to avoid large density gradients close 
to the origin. This simplifies our subsequent analysis and 
helps us get important contributions to the torque from 
just outside the gap, exactly where we want to test our 
theory. The factor 2.5 in the exponential is introduced 
for the parameter S to better correspond to the width of 
the inner, flat part of the gap. 

In the left panel of Figure Q] we plot the density profile 
for the disk model with parameters S — 5ft, E m ; n = 3 x 
10~ 4 , and with n = 2,4,6, respectively. These profiles 
represent rather broad gaps: for the models with S = bh 
shown in Figure Q] one finds that the density Soo/2 is 
reached at |x| k, 10h, \x\ « 8.5ft, and x « 7.7ft for 
n = 2, 4, and 6 respectively. Evidently, as n is increased 
the profile becomes boxier and the gradients at the gap 
edge (they are largest for 4ft < |x| < 15ft) grow. 

2.3. Shifted Lindblad resonances 

Non-zero pressure gradients intrinsic to inhomoge- 
neous disks affect the unperturbed azimuthal velocity ([?]) 
and change the locations of the Lindblad resonances, at 
which different harmonics of planetary potential couple 
to the disk, compared to the purely Keplerian case. In 
Appendix A we demonstrate that in the shearing sheet 
approximation a given k y harmonic has Lindblad reso- 
nance condition satisfied at xz(ky) satisfying the follow- 
ing implicit relation: 



ky(x L ) = ± 



2 k(xl) 



Xl 



ft 2 dlnE 



dx 



(24) 



where the modified value of the local epicyclic frequency 
k for the inhomogeneous disk is given by (see Appendix 
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blad resonance position given by Eq. i|24|] for the surface profile 
described by Eq. J23D with parameters S = 5h, S m i n = 3 X 10~ 4 , 
and n = 2,4.6. The dotted line indicates the uniform disk case 
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A) 



k(x) = ftp \ 1 + h 



,9 2 lnS c 
dx 2 



1/2 



(25) 



In a uniform disk (d In 'So/dx — 0) one recovers the usual 
expressions n — fl p and XL(k y ) — 2/(3k y ). 

For a class of density profiles (|23l) considered in this 
work the modification to the XL{k y ) relation (J24I) com- 
pared to the purely Keplerian case comes predominantly 
from the change of the local epicyclic frequency as given 
by equation (|25j) . For this reason we explore the depen- 
dence of k on x in some detail and plot it in the right 
panel of Figure Q] for the density profiles shown in the 
left panel of the same Figure. 

Deep inside the gap Eo is almost constant because of 
the assumed exponential dependence in equation (|23[) 
and k(x) — > Q p . As the gap edge is approached, at 
\x\ > Ah, <9 2 lnEo/<9x 2 increases and k does the same: 
in all displayed models k goes up to > 2Q, p . Beyond 
l x l ^ 5/i the density profiles change their curvature from 
positive to negative and k becomes smaller than f2 p . At 
even larger |x| the density profile flattens out and k con- 
verges to fip at |x| > 15ft.. The deviations of the epicyclic 
frequency compared to its value in a uniform disk be- 
come more extreme as we increase n (increase boxiness of 
Eo(x)), with the differences between models being more 
dramatic in regions where k < fl p . The minimum val- 
ues of k are 0.36fi p , 0.27fl p , and 0.03il p for profiles with 
n = 2, 4, and 6, respectively, see Figured! 

In fact, one can see from equation (f2"5) that depend- 
ing on the curvature of Eo the gap profile could have 
k 2 < at the gap edge, which would mean that the disk 
is Raylcigh unstable. For this reason in this work we ig- 
nore profiles with extremely deep gaps and strong density 
gradients at their edges since these are likely to result in 
k 2 < at some place. Our boxiest model with n = 6 
might be regarded as a limiting profile (for assumed val- 
ues of S and E m i n ) since the epicyclic frequency closely 
approaches zero at x ~ 6h. 
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Armed with this knowledge we plot in Figure [5] the 
dependence (|24| for the density model given by equation 
(|23|) with the same set of parameters as in Figure [T] 5 = 
5h, S m i n = 3 x 1CP 4 , and n = 2, 4, 6. By comparing these 
two Figures one can see that the resonances are displaced 
outwards relative to the uniform density case in regions 
where k > f2 p , while the opposite happens when k < 51 p . 
As a result, there is a concentration of the resonances at 
the gap edge where this transition occurs. 

Indeed, in a uniform disk with n — Q p the radial inter- 
val 4h < \x\ < Qh contains Lindblad resonances only for 
the modes with k y h in the range (0.11, 0.17). At the same 
time, for the gap in the form (|23p with n = 2 the same 
radial interval contains Lindblad resonances for modes 
satisfying 0.05 < k y h < 0.37; for n = 6 model modes 
with 0.009 < kyh < 0.4 are being excited in the same 
radial interval. Thus, concentration of resonances at the 
gap edge is more pronounced for density profiles with 
larger gradients meaning larger deviations of k from 51 p . 

Another striking feature of the XL{k y ) dependence 
clearly visible in Figure d is that even for not very sharp 
gaps (for rather low values of n) there are regions close 
to the gap edge where a single k y mode can be excited at 
two or three different locations (for a fixed sign of x), as 
happens e.g. for k y h = 0.3 in this Figure. This splitting 
of resonances in inhomogeneous disks has important im- 
plications for the torque behavior in Fourier space, see 

m 

3. NUMERICAL RESULTS 

3.1. Numerical procedure. 

Having specified the gap density profile the solution 
for the azimuthal velocity perturbation v is obtained by 
direct numerical integration of the equations (fTTj) . (|T2"j) . 
(116[) - (|20[) : knowing v we find u and (5E using equations 
d2U), (Hill- 
Numerical integration uses the same method as was 
employed in RP12 for uniform disks: we shoot nu- 
merical solutions from the origin and match them to 
the WKB outgoing waves fa r from the planet (see also 
iKorvcanskv fc Pollack! (J1993D ). As shown in NGG, the 
exact homogeneous solution to Equation (IT51) for a ho- 
mogenous disk, i.e. the parabolic cylinder functions, can 
be well approximated by the WKB outgoing waves 



v (x — > ±oo) 



^±i3k y x 2 /4h 



Of£yX 



(26) 



when x/h > (8/9) [l + (k y h) 2 ] /{k y h) 2 . Moreover, for 
this approximation to accurately represent an inhomoge- 
neous solution of Equation (fT5|) one requires xk y ^> 1, so 
that 4> < 1 and d<p/dx < 1. 

To additionally account for the presence of a density 
gap in the disk, we use the fact that for every profile 
dYi^/dx ~ » when x is several times the characteristic 
width of the gap 8 and, therefore, the asymptotic solution 
will still be given by (l2l)l) . All conditions together imply 
that x/h ^> mm{l/(k y h), 1,8/h}, which is the criterion 
we use to match our numerical solutions to the outgoing 
waves. 

Our numerical calculations use a 4-th order Runge- 
Kutta integrator with spatial resolution of h/AOO and 
x 6 [— 220h, 220h] for 720 uniformly log-spaced values of 



k y h between 0.01 and 15, and potential softening length 
of Wr 4 h. 

3.2. Density wake 

In Figure [3] we show one of the outcomes of our cal- 
culations — the behavior of the surface density pertur- 
bation Tii(x,y) in physical space. We obtain it by first 
calculating <5E using Equation (1221) with the numerically 
determined v and then performing the inverse Fourier 
transform. In Figure [3] azimuthal cuts of Ei(:c, y) at 
different values of x are normalized by (Eo(x)a;) 1 / 2 to 
ensure similar wake amplitude at different radial separa- 
tions from the perturber; this scaling is expected_| from 
the conservation of t he angular momentum flux carried 
by the density wave (jGoodman fc Ra fikov 200l|; iRafikovl 
I2002al) . 

In panel (a) we display our results for a shallow gap 
with a minimum surface density at the satellite position 
of O.lSoo (corresponding density profile is shown in panel 
(b)). For comparison we also display the wake profiles 
for a uniform density disk. Deep inside the gap, at x — h, 
where the Eo is essentially constant in our model we find 
our Ei to exactly coincide with the uniform disk calcula- 
tion, as expected (our normalization of Ei takes care of 
the difference in amplitude). However, further out from 
the planet the results start to differ. In particular, at 
x = lO/i the azimuthal location of the wake in a non- 
uniform disk gets shifted by « h relative to the uniform 
case, but the differences in the shape of the density per- 
turbation remain at the modest level. 

Situation is a bit different for deep gaps, as illustrated 
in panel (c) where we plot azimuthal density cuts through 
the wake for a gap profile with E m ; n = 0.0003Eoo (shown 
in panel (d)). Again, at x = h density profiles coincide 
for both uniform and non-uniform cases. However, at 
the gap edge (at x = 5h) there are now substantial dif- 
ferences between the profiles not only in the wake posi- 
tion but also in shape. The differences become more pro- 
nounced as the wave propagates even further: at x = lO/i 
the peak value of Ei is shifted by ~ Ah compared to the 
uniform case and the overall wake shape is considerably 
distorted. This is not surprising since the surface density 
eigenfunctions at the edge of the gap in a non-uniform 
disk are considerably distorted compared to their uni- 
form disk analogues, see i )5.1l 

Note that the amplitude of scaled Ei in Figure [3] does 
not vary too much with x. Since the angular momen- 
tum flux conservation predicts that Ei oc [Eo(a;)x] 1 ' 2 
it is clear that the amplitude of Ei grows as the wake 
propagates out from the bottom of the gap. At the same 
time, the relative density perturbation Ei/Eo(x) scales 
as [a;/Eo(a;)] 1 ' 2 and decreases as the wake climbs up the 
density contrast at the gap edge. This may have impor- 
tant implications for the wake damping, as discussed in 



4. ANGULAR MOMENTUM TRANSPORT: PRELIMINARIES 

One of the key characteristics of the disk-satellite cou- 
pling is the excitation torque density dT/dx (or equiva- 

7 This expectation is only approximate since the angular mo- 
mentum flux c arrie d by the wake can vary substantially across the 
gap edge, see i]6.2l Nevertheless, the proposed scaling works quite 
well in practice, which is obvious from Figure [3] 
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lently dT/dr) — the amount of the angular momentum 
added by the satellite tide to the density wave per unit 
radial distance x. Following RP12 we express it via the 
imaginary part of the density perturbation Q(<5S) as 

oo oo 

§-/*<-/(£).*■ <*> 

-oo y 

JTTi \ 

-^\ --4^fc y 03(^S]), (28) 

/ Ky 

where (dT/dx)k is the torque density contribution due 
to mode with wavenumber k y . Integrated (one-sided) 
torque is then defined as T(x) = J Q (dT/dx 1 ) dx' . Note 
that in the shearing sheet approximation there is no net 
torque and T(x) = —T(—x). 

Integrated torque is closely related to the angular mo- 
mentum flux (AMF) Fh{x) -- the amount of angular 
momentum carried by the wake, which can be written as 
(RP12) 



F H (x)= / F H ,k y {x)dk y , 



H.ku 



■■ 47rS (a;) [3fc(u)9fc(u) + S(v)3(u)] , (29) 



where Fn,k is the angular momentum flux carried by 
a particular azimuthal mode of the wave. Angular mo- 



mentum conservation demands the radial divergence of 
the AMF dF H ,k y /dx to be equal to (dT/dx) ky . RP12 
have demonstrated this to be the case for a uniform disk 
and their calculation can be trivially extended for an 
arbitrary non-uniform disk density profile. This is done 
(separately for each azimuthal mode) by writing Fu,k hi 
terms of v using (|2ip. differentiating with respect to x, 
and then rearranging the terms proportional to vdv/dx 
with the aid of equation (fT8"|) . As a result one indeed 
finds that dFn.k /dx = {dT / dx)k ■ 

4.1. Lindblad Resonance torque prescription 

To derive the asymptotic torque density behavior (U) in 
physical space GT80 have explicitly assumed that the an- 
gular momentum associated with a particular low-order 
(k y h < 1) harmonic of the perturbing potential gets 
added to the density wave exactly at the Lindblad reso- 
nance corre s pondin g to that mode, i.e xi,{k y ) — 2/(3fc a ). 
iDong et all (|2011al ) have extended this assumption to 
modes with arbitrary values of k y to obtain a so-called 
Lindblad Resonance (LR) prescription for the angular 
momentum flux Fjf^(x), which they used as a bench- 
mark for comparing with their numerical results for uni- 
form dis ks. This prescripti on is given by the following 
formula pong et al.ll2011af) : 



F h R (*) = l J d»j£^Fr B (») (*M 



.(*) 



FWKB (Ai) - 
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WKB 
H 



4 a (GM V 
3 \ c s 



x [2^(2/3) + ^ (2/3)] 2 , 



(31) 

where F H (n)/F^ KB (n) is the ratio of the actual AMF 
contribution carried by a particular harmonic with /i = 
k y h far from the perturber to the value of the same quan- 
tity obtained using WKB approximation. This ratio has 
been computed numerically by GT80 for arbitrary /i; in 
the limit /2 — >■ it tends to unity. 

In equation (|3"U1) n m i n {x) — (2/3)h/x is the mini- 
mum value of kyh for which the position of a corre- 
sponding Lindblad resonance (in a uniform disk) satis- 
fies XL(k y ) < x. Thus Fjf i (x) represents the sum of the 
AMFs (measured at x — > oo) carried by all individual 
Fourier modes having their Lindblad resonances lying at 
XL(ky) < x. Note that in computing fA m in(%) we do 
not use the positions of Lindblad resonances accounting 
for disk nonuniformity (Equation I24j) — we found this 
prescription for Xi,(ky) to provide Fjf 1 inconsistent with 
numerical results, see §6.11 

The prescription (|30|) is formulated for a disk with con- 
stant density Eoq. We extend it to disks with gaps by 
first defining a Lindblad resonance torque density in a 
uniform disk as dT hR (x)/dx\ u = dFjj R (x)/dx, and then 
using this torque density in equation ([3]): 



dT LR (x) _ E (aO dFJf{x) 



dx 



dx 



(32) 



Note that for |x| > h one has \Xmin ^5 1 an d 
dT JjR (x)/dx\ u reduces to the standard GT80 expression 
IJ3J). Most studies of tidal coupling in non- uniform disks 
have used this asymptotic form of the prescription (|32 ). 
i.e. equation ^ with dT/dr\ u given by equation ((5), 
to characterize the torque density; use of equation (|30p 
simply allows us to extend this prescription to arbitrary 
values of x. 

The prescription (|30[) is equivalent to assuming that 
each potential harmonic exerts torque only in the imme- 
diate vicinity of its Lindblad resonance, a notion that has 
been shown by RP12 to be incorrect in the shearing sheet 
approximation even for fj, < 1. In reality Lindblad reso- 
nances have finite width, which leads to the nontrivial in- 
terference between them, resulting in the negative torque 
density phenomenon at \x\ > 3.2h (Dong et al. 2011a; 
RP12) not captured by the GT80 analysis. Despite this 
failure of the assumed discreteness of the Lindblad reso- 
nance, we still use Fjf 1 ^) in this work to compare with 
our new results as it provides an interesting reference 
point. 

5. TORQUE BEHAVIOR IN FOURIER SPACE 

We start our investigation of the angular momentum 
exchange between the disk and the perturber by explor- 
ing the torque behavior for individual Fourier harmon- 
ics. We explore both the spatial structure of the torque 
density due to individual modes ( ij5.ll) . as well as their 
contribution to the full torque far from the perturber 
(» 

5.1. Spatial structure of(dT/dx) k 



In Figure U we show the torque density for individual 
Fourier harmonics (dT/dx)k (x) computed using equa- 
tion (|28p . at different radial separations x. These calcu- 
lations assume our fiducial density profile (|23|) with pa- 
rameters 5 = 5h, S m in = 3 x 10~ 4 , and n = 2. We show 
both the torque density per unit local surface density 
(i.e. (dT/dx)k y (x) normalized by Eo(x); left panels) as 
well as the pure (dT/dx)k (x) (divided by constant Eqo, 
right panels). We also display (dT/dx)k y \ u computed for 
a uniform disk (see Figure 3 of RP12), normalizing it by 
Soo in the left panels (which yields the torque density 
per unit local surface density in a homogeneous disk) 
and also multiplying by T,q(x)/'S^ >0 (to incorporate the 
local density correction) in the right panels. In all cases 
torque density is properly normalized to make it dimen- 
sionless. In the left panels of this Figure we indicate the 
position of the Lindblad resonances with arrows, both 
for an inhomogeneous (see equation [24j ) and a uniform 
disks. 

In panels (a) and (b), we plot (dT/dx)k v (x) for k y h — 
1. The corresponding Lindblad resonance lies at x/h — 
2/3, in the flat bottom of the gap and is thus the same 
for both uniform and fully non-uniform disks. Moreover, 
the waveform inside the gap perfectly matches that of the 
uniform disk since the profile has constant density there. 
Note that the torque density is not sharply localized at 
the Lindblad resonance position xr,(k y ) but is extended 
over a reasonably broad radial range around xl, with 
quite pronounced (even though damping) oscillations of 
(dT/dx)k {x) clearly visible even for x 3> xl,. This ad- 
ditionally confirms the finite and no n-negligible w i dth o f 
resonances previously pointed out bv lArtvmowiczl (|1993[ ) 
and RP12. 

At x > Ah, where Eo starts to increase, we observe that 
the torque density decreases its amplitude quite abruptly 
compared to the uniform disk calculation, see panel (b) . 
By looking at the shaded regions one can notice that the 
two calculations start differing at Eo ~ lCP 3 £oo where 
the density gradients become significant. This damp- 
ing of the amplitude compared to the homogeneous disk 
case is the consequence of incorporating the disk non- 
uniformity in the fluid equations, also clearly visible for 
other values of k y h at the gap edge (see below). 

In panels (c) and (d) we show the torque density for 
kyh = 0.32 mode, for which there are multiple Lind- 
blad resonances at x/h — 2.1, 3.8, and 4.5, see Figure 
[2] (the uniform disk has, of course, just a single one at 
x/h = 2/(3k y h) « 2.1). These are indicated with arrows 
in panel (c) and it can be seen from both panels that 
the extra resonances located at x ~ Ah give rise to pos- 
itive torque density contribution and overall significant 
phase shift compared to the uniform disk case. Also, the 
amplitude of the torque density oscillations far from the 
resonances is greatly reduced in the self-consistent cal- 
culation of (dT/dx)k (x) so that in panel (d) we even 
multiply the uniform disk curve by 0.05 for it to have an 
amplitude comparable to that of the fully non-uniform 
disk calculation. 

The mode with k y h « 0.14 shown in Figure |4^,f has a 
single Lindblad resonance at xl ~ 5h, which coincides 
with the resonance position in a uniform disk (see Figure 
[2|). Despite this coincidence the waveforms for the uni- 
form and non-uniform disk calculations are very different 
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0.001 x Eoo, 0.01 x Eoo, 0.1 x Eoo, and 0.5 x Eoo. 



from each other both in terms of shape and the overall 
amplitude (note that the uniform disk curve in panel (f ) 
has been multiplied by 0.1). 

The mode with k y h — 0.08 shown in Figure 2g, h has 
its (single) Lindblad resonance at x « 5.1h, considerably 
shifted inward from the corresponding uniform disk res- 
onance position at x = 8.3ft. The torque density normal- 
ized by Eo(x) (see panel (g)) has a rather broad spatial 
distribution extending from x ~ 5ft so x ~ 8ft before 
starting to oscillate at an amplitude reduced roughly by 
a factor of 2 compared to the uniform disk calculation. 

Finally, in panels (i) and (j) we show the mode k y h = 
0.04 for which xj, ~ 16ft. lies where the density gradients 
are negligible. As a result, calculations in a uniform and 
nonuniform disk cases perfectly coincide, demonstrating 
that far from gap torque excitation is well described by 
the uniform disk theory of RP12. 

To summarize, the spatial distribution of the torque 
density for individual Fourier harmonics can be substan- 
tially modified in non-uniform disks in regions where the 
density gradients are large, compared to the case of a 
constant density disk. The differences arise both because 
of the displacement of the Lindblad resonances and be- 



cause of the different mathematical structure of the fluid 
equations in the non-uniform disk case. 

5.2. Angular momentum flux in Fourier space 

We now study the behavior of the AMF as a function 
of k y for x — > oo. Since, as mentioned in [JH FH,k y (x) 
is identical (up to a constant offset) to the integrated 
torque we look at the behavior of 



™«'=/l§ 



fix 



(33) 



in Fourier space, where (dT '/ ' dx)k is given by our numer- 
ical calculations. 

As a point of comparison we also use the variation of 
Fniky) = Fn,k y {x — > oo) with k y in a homogeneous 
disk (previously computed by GT80) which we addition- 
ally multiply by Yio(xL,{ky)) to account for the disk non- 
uniformity (where we take XL{k y ) — 2/(3k y )). This re- 
scaling implies that each mode carries the AMF from a 
uniform disk corrected by the local density at the posi- 
tion where the excitation of this mode occurs in a uni- 
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Fig. 5. — Angular momentum flux (or integrated torque) far from 
the planet (x — ¥ oo) in Fourier space (logarithmic coordinate s at 
the top, linear at the bottom) for a gap model given by Eq. (123 I t 
with <5 = 5h, E m ; n = 3 X 10~ 4 , and n = 2. The black solid line 
indicates the results from equation (133 I t and the red dashed line 
results from the LR theory described in JJ5.2I 



form disk. This procedure is a version of the LR theory 
described in i M.ll in Fourier space, as it incorporates the 
two major ingredients from that prescription: the local- 
ized mode excitations at a Lindblad resonance and the 
AMF computed for a uniform density disk. 

In Figure [5] we plot the torque from our calculations 
and the AMF from the LR theory for our fiducial model 
with S — 5h, S m i n =3x 1CT 4 , and n = 2. The first thing 
we notice is that both expressions coincide in the limits 
of long wavelengths k y h < 0.05 and short wavelengths 
k y h > 2. This is expected to happen since in both limits 
the satellite excites waves in the flat parts of the disk — 
far outside the gap for k y h < 0.05 and deep inside the 
gap, in the flat part of the density profile for k y h > 2. 
For that reason the positions of the resonances and the 
waveforms coincide in our non-uniform calculation and 
in a uniform disk. But in the intermediate range of k y 
there are appreciable differences between the LR theory 
and our calculations, mainly due to effect of the shifted 
Lindblad resonances as we discuss next. 

First, for k y h ~ 0.05 — 0.1 Figure [2] demonstrates that 
Lindblad resonances are shifted towards the perturber 
and the excitation occurs in lower density region of the 
disk compared to the case when resonances lie at x = 
2/(3k y ) (as appropriate for homogeneous disk). Because 



of this reduction of the background surface density the 
torque produced by modes in this range of k y ends up 
being lower than the AMF computed according to the 
LR theory (which assumes xl = 2/(3fc a )), even though 
the coupling to the perturbation potential is stronger at 
smaller x|. Conversely, modes with k y h ~ 0.1 — 0.2 have 
their resonances shifted outwards, into the higher density 
part of the gap profile, which explains their higher values 
of T{k y ) compared to the LR theory AMF. 

At higher values of k y there is a point where reso- 
nances split and take place at multiple locations, as dis- 
cussed in §2.31 According to Figure [2] this is the case for 
k y h > 0.2, which explains the intricate shape of T(k y ) 
for these modes — a peak at k y h = 0.2 — 0.3 where the 
Lindblad resonances are split in three (see Figure Sfc,d 
for illustration of the torque density behavior for these 
modes). The innermost resonance almost coincides with 
that in a uniform disk, but the other two are located fur- 
ther out, in a higher density region of the gap, explaining 
the spike of T(k y ) for these values of k y . 

Complicated structure of the T(k y ) curve around the 
peak at k y h ~ 0.32 does not have a straightforward ex- 
planation and is probably related to the intricate behav- 
ior of the torque density when the excitation is taking 
place at multiple locations at the gap edge. 

6. TORQUE BEHAVIOR IN PHYSICAL SPACE 

We now explore the behavior of the torque integrated 
over all Fourier modes in physical space. We look at the 
spatial distribution of the torque density in i j6.ll and of 
the integrated torque T(x) in i|6.2l 

6.1. Torque density 

In the top panel of FigureEJ we plot the torque density 
dT/dx (integrated over all k y , see equation [25] ) per unit 
of the local surface density E (x) for our fiducial model 
and make a comparison with the LR theory behavior 
given by equation (1321) . In the lower panel of Figure [6] we 
directly compare dT/dx and dT hB -/dx. One can clearly 
see a significant discrepancy between the two. 

The disagreement at \x\ < 5h, where the surface den- 
sity is close to constant (the bottom part of the gap) is 
expected based on the RP12 results — one can see that 
our dT/dx is negative for x ~ (3 — A)h which is just the 
manifestation of th e negative torque density phenomenon 
(jDong et al.M2011af > not captured by the dT LR /dx. Also, 
the overall shape of dT/dx and dT hR /dx is different for 
|x| < 3h — again well kno wn from the uniform disk stud- 
ies of iDong et all (|2011al ) and RP12. 

The discrepancies at \x\ > 5h where the density gra- 
dients are important can only be understood by fully 
accounting for the disk non-uniformity. According to 
GT80 for |x| > h the LR theory predicts H^dT/dx « 
2.5(/i/x) 4 (in this section we refer to the torque den- 
sity in units of Y, O0 G 2 M 2 /{hc s ) 2 ), see equation ((H) and 
i j4.ll However, our self-consistent nonuniform disk cal- 
culations suggest that the torque density decays expo- 
nentially oc exp (— 0.8x/h) for |x| > 5h, as shown in the 
upper panel of Figure HO This result casts serious doubt 
on all torque prescriptions used in the literature that 
adopt the GT80 functional form to account for the tidal 
torque. The exponential fall-off is clearly related to the 
presence of strong density gradients at the gap edge but 
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Fig. 6. — Torque density dT/dx obtained self-consistently for 
a non-uniform disk using our nume rica l results (solid black lines) 
and that from the LR prescription f £|4.1l dashed red lines) for a gap 
model given by Eq. I I23I) with parameters S = 5h, S m ; n = 3 X 10 — 4 , 
and n = 2. In the upper panel we normalize by the surface den- 
sity T,o(x) and add a numerical fit to show the exponential fall-off 
of dT/dx outside the gap. In the lower panel we normalize by 
the density at infinity Soo and also display dT/dx computed using 
the prescription in Eq. 10 with the uniform disk torque density 
properly accounting for the negative excitation torque density phe- 
nomenon (RP12). The latter curve demonstrates the conceptual 
failure of a simple prescription l T3t . 

we could not find a simple qualitative explanation for this 
particular form of the torque decay. We just note here 
that this behavior of dT/dx has nothing to do with the 
presence of the exponential term in equation (j2"3")) , which 
is introduced only to guarantee the flat density profile in 
the bottom part of the gap. We provide more details on 
the exponential decay of dT/dx in Appendix |E] 

Apart from the discrepancy in the functional form, the 
self-consistent specific torque density is also more con- 
centrated towards the perturber at the gap edge, and 
has a larger amplitude there. This effect can be under- 
stood from our previous analysis of the shifted Lindblad 
resonances that tend to accumulate at the gap edge (see 
Figures [2] and 0|) , increasing excitation torque density in 
this region and reducing it outside of it, at |x| > 10ft. 
This effect has not been taken into account in previous 
studies but it does drastically change the overall shape 
of the torque density. 



We tried to improve the performance of the LR the- 
ory fi )4.1[) by assigning torque contributions of individual 
modes not to xl — 2/(3^) (which is appropriate only for 
a uniform disk), but to XL\k y ) given by the self-consistent 
calculation in a non- uniform disk, see equation (|24j) . Un- 
fortunately, the agreement with our fully self-consistent 
calculation has gotten even worse (for this reason we did 
not pursue this idea further), which suggests that prop- 
erly accounting for the resonance overlap is also very im- 
portant for getting the correct spatial distribution of the 
torque. 

Far enough from the gap edge, where So (a;) flat- 
tens out, the excitation torque density becomes neg- 
ative, which is a direct m anifestation of the negative 
torque density phenomenon (D ong et al.ll20lia() . For the 
adopted gap profile this happens at x w 15ft., and dT/dx 
stays negative out to infinity and agrees with the asymp- 
totic behavior — 0.63(/i/x) 4 derived in RP12. 

Additionally, in the bottom panel of Figure [5] (dot- 
dashed curve) we plot dT/dx calculated according to the 
simplified prescription ([3]) but using the correct excita- 
tion torque density for the uniform disk dT/dx\ u calcu- 
lated in RP12 instead of the GT80 prescription (j4|). One 
can see that this calculation fails miserably compared to 
the correct dT/dx: it predicts negative excitation torque 
density for x > 3. 2 ft,, which is the direct consequence of 
the negative torque density phenomenon in the uniform 
disks. Because of this, as we will later see in Figure[7l the 
use of this prescription results in the negative integrated 
torque deposited in the density wave by the perturber, 
which cannot be true in a real disk. 

The point of this last exercise was to illustrate the 
failure of the simple prescription Q, which yields un- 
physical results even when it uses physically motivated 
ingredients — the spatial behavior of dT/dx\ u derived in 
RP12 as opposed to the asymptotic formula (j4|). 

6.2. Integrated torque and AMF 

In Figure [3 we plot the integrated torque T(x) for 
our fiducial model as well as the AMF Fu{x) calculated 
using equation (f2"5)l . The two have small relative vertical 
offset but are otherwise identical, as expected from the 
discussion in fj4j We also display T LR (x) obtained by 
integrating dT LR /dx given by equation ((32|) from zero to 
x. 

All these expressions almost coincide within the gap. 
This happens because in the flat-bottomed gap different 
methods of computing integrated torque should yield the 
same result: a total torque (or AMF) of w 0.93 (GT80) 
in units of the plot times Emm = 3x 10~ 4 . This equiva- 
lence of the full torque calculations in the Fourier space 
(T LR (cc)) and physical space (T(x)) has been previously 
mentioned in RP12. 

At the gap edge, for |a;| > 5h, one finds that T(x) in- 
creases faster than T LR (a;), which is a direct consequence 
of the torque density dT/dx in a self-consistent non- 
uniform disk calculation being rather concentrated to- 
wards the perturber, see H6.ll Quite remarkably, despite 
this difference at the gap edge the accumulated torque 
at infinity nearly coincides with that given by the simple 
theoretical prescription T LR (x): our self-consistent cal- 
culation yields full torque w 1.2 x 10~ 3 , while the the LR 
theory predicts ~ 1.3 x 10~ 3 (see Figure[7]at x — 30ft. for 
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Fig. 7. — Integrated torque T(x) from our numerical results (solid 
black line), the LR theory (dashed red line), and the integration 
of the torque density from uniform density disk (RP12) using the 
prescription $3% (dot-dashed blue line). Das hed green line displays 
the AMF Fjj{x) obtained from Equation J29H and, as expected, 
agrees with T(x). Note that the RP12 prescription for dT/dx\ u 
used in combination with equation J3} resul ts in negative integrated 
torque as x — > oo, which is unphysical (sec i|6.1l for details). Results 
are for a gap model i'2'M ) with 5 = 5h, S m ; n = 3 X I0~ 4 , and n = 2. 

reference). This means that the total torque exerted on 
the disk (but not the spatial distribution of the torque 
density!) is well described by the simple theoretical pre- 
scription (essentially T hR (x) based on asymptotic scaling 
[4]) that has been broadly used in the literature. 

We do not have a fully satisfactory explanation for this 
coincidence except for the following observation. In the 
lower panel of Figure [5] we display the torque in Fourier 
space (on linear scale) for k y h < 0.4. There it can be seen 
that our non-uniform disk calculation generate almost 
the same amount of torque as the theoretical prescrip- 
tion (integrals under the curves are roughly the same), 
meaning that the torque deficit with respect to the LR 
theory visible for k y h < 0.1 is almost exactly compen- 
sated by the torque excess for k y h > 0.1. This means 
that the effect of resonances, which are being shifted to 
lower density regions resulting in lower torque contri- 
bution (smaller k y h) is closely counterbalanced by the 
torque enhancement due to the resonances displaced to 
higher density region (larger k y h). This is an interest- 
ing result, which only weakly depends on the exact gap 
profile as we empirically show next. 

7. EFFECT OF VARYING GAP SHAPE PARAMETERS 

So far, we showed the results obtained for a specific gap 
density profile given by equation f|23[) with a certain set 
of parameters, namely S — 5ft., S m in = 3 x 10~ 4 , and n = 
2. It would of course be much better to discuss torque 
structure not for an arbitrary model of the gap profile 
but for the one which is realized in nature. However, the 
calculation of such profile must self-consistently couple 
the calculation of the torque excitation (such as the one 
done here) with the description of the wave damping, and 
this is beyond the scope of this work. However, we can 



still explore the general trends in the torque behavior as 
the various properties of the gap (its depth, width, etc.) 
are varied. This is what we do now. 

7.1. Steepness of the gap edge (variation of n) 

As discussed in S J2.2I and shown in Figure [TJ the pa- 
rameter n controls the steepness of the density gradi- 
ents at the gap edge. In panels (a) and (d) of Figure [3J 
we plot the excitation torque density and the integrated 
torque, respectively, for n = 2,4,6, while fixing 5 = 5h 
and £ min = 3 x 10" 4 . 

As n is increased the gradients become larger and, 
therefore, the amplitude of the variation of k 2 around 
the gap edge becomes larger. This leads to higher reso- 
nance accumulation closer to the perturber, and So gets 
higher there as well. As a result, the torque density gets 
more concentrated towards the gap center with increas- 
ing n, see Figure [5^: for n = 6 almost all the torque is 
excited within the range 5ft < x < 10ft, while for n = 2 
this region extends out to x ~ 15ft. 

The integrated torque also increases with n since closer 
to the gap center excitation by the perturber is more 
effective, depositing more angular momentum into the 
density waves. We see from Figure [SJi that the density 
profile with n = 2 results in a total torque of ~ 1.2 X 10 -3 
in units of (GM p )' 2 'E, O0 /hc 2 ., while using n = 6 we obtain 
-2.1xl0- 3 . 

This Figure also demonstrates that for all values of n 
the LR theory predicts more extended torque density dis- 
tribution than we find in self-consistent calculations, see 
£16.11 Nevertheless, the integrated torque at \x\ —> oo is 
still well predicted by the LR theory, the point we pre- 
viously made in tj6.2l The largest relative difference of 
~ 13% is found for n = 6, for which we obtain the inte- 
grated torque of ~ 2.1 x 10~ 3 , while the LR prescription 
gives — 2.4 x 10~ 3 . 

7.2. Width of the gap (variation of 5) 

In Figure [5b, e we vary the width of the gap by slightly 
changing 5, but fixing £ m j n = 3 x 10~ 4 and n = 2. As ex- 
pected, for narrower gaps the torque density peaks closer 
to the planet and the integrated torque reaches higher 
values, simply because the strength of the tidal coupling 
is higher closer to the perturber. 

Figure [SJ3 shows that the integrated torque predicted 
by the LR theory tends to better match the self- 
consistent calculation as the gap width is increased. Also, 
the dependence of the T(x — > 00) on S is not linear: de- 
creasing S from 6/1 to 5.5/i results in ~ 14% increase of 
the integrated torque, while going from 5.5ft to 5ft re- 
sults in the increase of — 25%. Using the fact that the 
LR theory gives a reasonably good approximation to the 
integrated torque behavior, and adopting the prescrip- 
tion (|3|) with dT/dr\ u given by equation (|4]) one can es- 
timate that the total torque excited outside the density 
gap scales as T(x > 6) oc (8/h) 3 . We verified this simple 
scaling by calculating integrated torque for density pro- 
files with S = (4 — 8)ft. Thi s behavior is in agreement 
with lPapaloizou fc Lin! (|1984D who present the same scal- 
ing for the one-sided torque. 

7.3. Gap depth (variation of S m j n/ ) 
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Fig. 8. — Torque density (upper panels a, c, and e) and integrated torque (lower panels b, d, and f) for a gap model given by Eq. I|23ll 
varying the three parameters n, <5, and £ m i n one at a time with respect to our fiducial choice 6 = 5/i, E m i n = 3 X 10 -4 , and n = 2. In 
panel (a) and (d) we vary n, fixi ng <5 , S m i n , in panel (b) and (e) we vary S, and in panel (c) and (f) we vary S m i n . Dashed line indicates 
the results from the LR theory f i|4.U , 



In Figure [8fc,f we vary the dimensionless gap depth 
Smin from 10 -4 to 10 -3 , while fixing n — 2 and 6 — 5h. 
One can see that changing £ m i n only affects the torque 
density and the full torque accumulated inside the gap 
(in a linear fashion) ; outside the gap (for x > 4/i) torque 
structure remains unchanged. 

As Smin is increased there is a point at which the 
integrated torque starts being dominated by the exci- 
tation inside the gap, in the immediate vicinity of the 
perturber. One can easily predict when this happens. 
For a flat-bottomed gap profile like the one we consider 
in this work the torque accumulated inside the gap is 
T(x < 6) ~ 0.9E min 'E oo G 2 M$/hc 2 s (RP12 and GT80) if 
the gap width d is larger than about 2h — the extent or 
the region where most of the torque excitation occurs in 
a uniform disk. This exceeds the torque T(x > 5) excited 
outside the gap for 



> 



T(x > 8)hc 2 s 



0.9£ c 



D G 2 M p 2 



2.8/i J 



E (x) 



dx, (34) 



since the value of T(x > 5) is well approximated by the 
LR theory using the prescription ([3]) with dT/dr\ u given 
by equation (j4|). In the case shown in Figure [8b, f the 
critical value of E m i n is rs 10~ 2 , i.e. the integrated torque 
in this case is dominated by the excitation at the edge of 



the disk only when the density inside the gap is depleted 
below 1% of Eoo. 

Finally, after trying different combination of parame- 
ters for the density profile, we observe that the torque 
density always decays exponentially outside the gap in 
regions where the density gradients are important, as 
discussed in ij6.ll In Appendix B we construct an empir- 
ical fit for the torque density in this region represented 
by equation (|B1|) . 

8. DISCUSSION 

The main result of this work is a new, self-consistent 
calculation of the torque excitation by a massive per- 
turber in a non-uniform disk. We clearl y show that the 
existing torque density prescriptions (ILin fcP apaloizou 
1986t iTrilling et all 119981: lArmitage fe Natarajanl 12002 ' 



02 



lArmitage et al.l l2002| ) . majority of which are based on 
the GT80 asymptotic result (|3]), provide inadequate de- 
scription of the spatial distribution of dT/dr in disks with 
gaps. 

In particular, we find that for deep gaps the torque 
excitation is more concentrated towards the perturber 
than the standard prescription ([3]) with dT/dr\ u from 
equation (UJ) would predict. Also, in regions where den- 
sity gradients are non-zero we find the fall-off of the 
torque density to be exponential as opposed to the a; -4 
dependence that is assumed in all previous works (see 
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Papaloizou fe Lin! (J1984J ); ILin fe Papaloizoul (|1986l ) and 
references therein). Moreover, far from the gap edge, 
where the density flattens out, we recover the nega- 
tive torque densit y phenomenon previously found by 
iDong et al.l (|2011a[ ) and RP12 in a uniform disk. 
There are three main reasons for these differences: 

1. Lindblad resonances in non- uniform disks are 
shifted with respect to their constant density disk 
analogues, and often occur at multiple locations 
(» 

2. Spatial behavior of the eigenmodes of the per- 
turbed fluid variables is different in a non-uniform 
disk (g53J). 



3. Interference of different Lindblad resonances im- 
portant for both uniform and non-uniform disks is 
usually not accounted for by the LR-type prescrip- 
tions ( gup . 

One can try to improve the existing semi-analytical 
torque prescriptions. Previously there have been some 
attempts to incorporate the (weak) disk non-uniformity 
in the LR prescriptio n for dT/dr by accounting for 
only t he resonance shift ( Menou fe Go odman 2004; iWardl 



119971 iMatsumura et al.l I2007F . We checked the perfor 
mance of this procedure in £ )6.1I but found that it does 
not improve (at least for strong density variations) the 
agreement with the self-consistent calculation. And the 
other two issues listed above arc clearly very difficult to 
account for in any simple way. 

Unfortunately, an accurate calculation of dT/dr for an 
arbitrary So( r ) profile, which is necessary for studies of 
the couple d evolution of the p crturber and the surr ound- 
ing disk (ILin fc Papaloizoul 119861: iTrilling et alj 119981: 



Ivanov et all 119991: lArmitage et al.l I2002t IChang et al.l 
2010), is rather computationally-intensive as outlined in 



N3- 11 In this regard it would be nice to have a simple ana- 
lytical fit to the actual excitation torque density behavior 
and in Appendix [B] we provide a simple one-parameter 
formula (|B1[) for dT/dr for exactly this purpose. Even 
though there is some ambiguity in choosing the proper 
fit (it depends on a single parameter a, which is some- 
what different for different gap profiles) this approach 
provides a basis for finding better semi-analytical rep- 
resentations of the dT/dr behavior useful for long-term 
numerical calculations of the disk-satellite interaction. 

Despite the different predictions for the torque density 
dT/dr we found that the LR theory outlined in fc|2.3l does 
produce a reasonable description of the total integrated 
torque, within ~ 20% of that found in a self-consistent 
calculation for diff erent values of the gap d ensity profile, 
see SJ7] Previouslv lPapaloizou fc Linl ( 19841 ) carried out a 
numerical torque calculation for a polytropic disk trun- 
cated by a massive perturber located at the distance 5 
from the disk edge, and found that the integrated torque 
scales as ~ <5 -3 . This is in agreement with our results 
described in J J7.2I and supports the observation that the 
LR theory provides reasonably accurate description of 
the full one-sided torque exerted on the disk. 

Based on this remarkable result one may decide that 
the orbital evolution of the perturber should be insen- 
sitive to the discrepancy in the torque density calcula- 
tion described above, simply because the variation of the 



perturber's angular momentum depends on the full inte- 
grated torque it exerts on the disk. This, however, is not 
true since the integrated torque does sensitively depend 
on the gap density profile (see [J7]), and the latter is a 
function of the torque density distribution in the disk. 
Thus, dT/dr is in fact involved (even though indirectly) 
in determining the rate of migration of the perturber. 
This makes any statements regarding the Type II migra- 
tion speed in protoplanetary disks and the rate of SMBH 
inspiral due to interaction with the circumbinary disk 
dependent on our understanding of the torque density in 
non-uniform disks, which our work aims to provide. 

8.1. Nonlinearity of the density wave 

We noticed in £|3.2l that the conservation of the angular 
momentum forces the relative amplitude of the density 
wake Si/Eo to decrease as the wave climbs up the density 
gradient at the gap edge. Interestingly, even though the 
normalized Ei/Eo does not vary much as the wave travels 
through the gap edge, the angular momentum density 
varies a lot; this is because the latter is sensitive to the 
overall wake shape and not just its amplitude. 

The variation of the relative density perturbation has 
interesting implications for the strength of the density 
waves outside the gap. In particular, from the fact that 
the normalized wake amplitude in Figure [3] does not vary 
much with x it follows that the wave nonlinearity mea- 

1 /2 

suredby the ratio Si /Eo oc (x/T,q(x)) should decrease 

by a factor of ~ [(<5//i)/E m i n ] between the bottom 
of the gap and the part of the disk just outside of the 
gap, where Eo ~ Eoo. For clean gaps with small S m i n 
this factor can be quite significant, « 0.2 for a gap with 
the width of S = 5h and density contrast at the bottom 
E m in = 10~ 2 (this factor is ~ 0.05 for the situation shown 
in Figure [3];). Thus, even if the density wave starts out 
highly nonlinear (Ei/So > 1) inside the gap near the 
planet, it may become linear upon propagating out of 
the gap. The initial wave nonlinearity in the flat part of 
the density profile near the planet is Ei/Eo ~ M p /M t h, 
where M t h = c|?/(S!G) is the characteristic planetary 
mass at which the Hill radius becomes comparable to 
the disk scale height. This mass is about IOM0 at 1 AU 
and about 6OM0 at 10 AU. Thus, a Jupiter mass planet 
at 10 AU would generate a highly nonlinear wave with 
S1/E0 ~ 5 at the bottom of the gap but the nonlinearity 
should drop to £i/Eo ~ 1 outside the gap for S = 5h and 
E m i n = 10~ 2 , even if we neglect the wave damping (not 
captured by our linear theory). This demonstrates that 
even rather massive planets capable of opening gaps may 
still have their density waves in the linear regime outside 
the gap, depending mainly on the density contrast S m j n . 
This point is rather important for the possibility of the 
nonlinear wave damping. IGoodman fc Rafikovl (2001) 
have studied the nonlinear density wave evolution in a 
uniform disk and demonstrated that shock formation re- 
sulting from nonlinear effects can lead to effectiv e damp- 
ing of t he wave. This study has been extended by|Rafikov 
(2002a) to include the possibility of radial variation of 
Eo- Our linear results suggest, in agreement with this 
latter work, that the nonlinear wake distortion (which 
ultimately results in shock formation) should slow down 
as the wave gets out of the gap, shock formation gets 
postponed and the nonlinear dissipation becomes less ef- 
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fective at transferring wave angular momentum to the 
disk material. 

For massive planets the wave might shock close 
to the planet, while it is still propagating through 
the roughly uniform bott om part of the gap profile: 
IGoodman k Rafikov (2001) estimate the shocking length 
of just 2ft, for a 2M ffi planet at 1 AU, and it is even smaller 
for more massive planets. Then the efficient density wave 
dissipation will start inside the gap but then will appre- 
ciably slow down as the wave climbs out of the gap (even 
though the shock will still persist). This complication 
should certainly affect the picture of the nonlinear wave 
damping and may produce important feedback for the 
self-consistent calculations of the gap density profile. 

8.2. Astrophysical implications 

Theory of disk-satellite interaction has been exten- 
sively used for understanding the interaction of proto- 
planetary disks with embedded planets and the evolu- 
tion of SMBH binaries surrounded by circumbinary ac- 
cretion disks, among other astrophysical settings. Our 
extensions of this theory have direct impact on our un- 
derstanding of these systems. 

In particular, our finding of the torque density being 
quite concentrated near the perturber suggests that re- 
gardless of the details of the wave damping, the angu- 
lar momentum carried by the density waves should be 
deposited in the disk closer to the perturber than the 
conventional theory would suggest. This certainly facil- 
itates the formation of density gaps or cavities since for 
the same mass of the perturber a narrower annulus of 
the disk needs to be cleaned by the tide. This may make 
Type II migration possible for lower mass planets than 
was previously thought. 

To fully understand the impact of our results on the 
Type II migration and SMBH inspiral one needs to know 
a self- consistent gap shape, which can be computed only 
when the wave damping mechanism is specified. In this 
regard we note that the majority of studies of the gap or 
cavity opening simply ignore the issu e of wave damping 
process and take dT/dr \d = dT/dr ()Lin k Papaloizoul 
ll986tlChang et aill2010( ). i.e. assume immediate deposi- 
tion of the excitation torque into the disk. However, the 
validity of this approximation has never been demon- 
strated, and calculatio ns that do properly take w ave 
damping into account (jWardl 119971 : iRafikovl l2002bH re- 
sult in a rather different picture of the gap opening. Our 
results on the wave amplitude evolution and its implica- 
tions for the nonlinear wave damping in ^3.21 18.11 should 
be useful for understanding the spatial distribution of 
dT/dr\d and calculating the gap shape (and orbital evo- 
lution of the perturber) fully self-consistently. 

Our present work is cast in terms of the shear- 
ing sheet approximation, which limits the direct ap- 
plication of its results to systems in which gaps are 
narrow compared to the semi-major axis of the per- 
turber. This i s the case for e.g. planets in proto- 
plane tary disks (|Lin k Papaloizoulll986t iTakeuchi et al.l 
1996) and SMBH binaries with extreme mass ratio sur- 
rounded by an accretion disk (so-called extreme mass- 
ratio inspirals (EMRIsl) (lArmltage &: NataraiaiJ l2002t 
IChang et al1l2010t iKocsis et alJl2011li ~ 

SMBHs binaries of similar mass do not fall in 
this category since they are known to clear out a 



large and clean cavity in circumbinary gaseous disk, 
with inner radius comparable to t he binary separation 
(e.g. J MacFad ven k Milosavlievicl (120081). ICuadra et all 
(20091 iShi et all (|2011h . IFarris et al.l (|2011h ). However, 
even though the shearing sheet approximation breaks 
down in these systems, our results can still be used to 
gain qualitative insight into the dynamics of cavity clear- 
ing. In particular, concentration of the torque density to- 
wards the perturber means that cavity should be opened 
by lower mass ratio SMBH than was thought before. 

Also, our results can help understand the excitation 
torque density patterns derived from simulations, most of 
which (IMacFadyen k Milosavlievicl I2008t ICuadra et aT 



2001 iShi et all I201H IFarris etafl [2oTlt iRoedig et al 



2012f ) agree that outside the cavity dT/dr exhibits com- 



plicated oscillatory pattern (some of these studies include 
additional physics such as MHD or general relativity). 
This behavi or was not recognized be f ore an d some earlier 
studies (e.g. lArmitage k Nataraianl (|2002l )) tried model- 
ing numerically derived dT/dr in the form similar to the 
GT80's result (@J) but with the amplitude reduced by a 
factor of ~ 10~ 2 . Such low amplitude most likely results 
from averaging actual spatially rapidly oscillating pat- 



tern of dT/dr and then applying the ex 



scaling 



to reproduce the total integrated torque in the simula- 
tions. Clearly, usage o f this prescription is likely to lead 
to m isleading results (jliu k Shapiro! 120101 IChang et all 
2010) in modeling gap shape. 

IMacFadyen k Milosavlievicl (|2008ft suggested that the 
oscillatory behavior of the torque distribution in their 
simulations is consistent with forcing at the 2:3 (m = 2) 
outer lindb lad resonance and tried fi t ting linear wave- 
forms from iMever-Vernet k Sicardvl (J1987I ) to spatial 
variation of the perturbed fluid variables. Even though 
they were able to reproduce qualitatively the oscilla- 
tory behavior of the numerical dT/dr the amplitude and 
phase of the torque density were substantially different 
between the numerical calculation and the analytical wa- 
verform, qualitatively consistent with the lindblad res- 
onance shift and reduction of perturbation amplitude in 
our calculations, see Figure 4f. This is certainly not sur- 
prisin g since not only were the IMever-Vernet k Sicardvl 
(|1987h waveforms derived for a uniform disk, they are 
also valid only in the close vicinity of the lindblad reso- 
nance. The latter limitation was removed in RP12 who 
came up with a fully global analytical approximation for 
the waveforms, but again only for a uniform disk case. 
Our present work is free from these constraints since 
it provides a way of computing (dT/dr)k for arbitrary 
disk surface density profile at any distance from the res- 
onance, see §5.11 Future extensions of our approach to 
full cylindrical geometry will be directly applicable to the 
astrophysical systems such as SMBH binaries in which 
broad gaps or cavities are typical. 

9. SUMMARY 

In this work we have carried out a linear study of the 
tidal disk-satellite interaction in the shearing sheet ap- 
proximation assuming the disk to have a non-uniform 
density structure to model the effect of density gaps 
around the perturber. The disk non-uniformity is self- 
consistently incorporated in the fluid equations, includ- 
ing the modification of the disk rotation curve due to 
pressure gradients, which are especially prominent at the 
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gap edges. This allows accurate calculation of the pri- 
mary fluid variables such as the density perturbation and 
torque distribution in both physical and Fourier space. 

Because of the rotation curve modification the Lind- 
blad resonances get shifted towards the gap edge, often 
occurring at two or three locations on one side of the 
perturber. Because of this and the terms in fluid equa- 
tions related to density gradients the waveforms of the 
modes excitated close to the gap edge are considerably 
modified compared to those obtained in a uniform den- 
sity disk. We find the torque density in physical space 
to be more concentrated towards the perturber and dif- 
ferent from the uniform disk theory predictions typically 
by a factor of several in this region. At least for the par- 
ticular gap profiles considered in this work we find that 
the torque density at the gap edge drops exponentially 
with the distance from the perturber (as opposed to the 
power law scaling expected from uniform disk theory). 
In parts of the disk where the density is roughly uniform 
we observe the negative torqu e density ph e nomen on in 
agreement with the results of iDong et al.l (|2011al ) and 
RP12. Despite these differences in the torque density 
distribution the total angular momentum flux driven by 
the perturber through the disk seems to agree with the 



existing prescriptions at the level of ~ 10%. The relative 
perturbation of the surface density is shown to go down 
in amplitude as the density wave propagates out of the 
gap, which has important implications for the nonlinear 
wave evolution and damping. 

These results suggest that the process of gap opening 
in protoplanetary disks and gas clearing around SMBH 
binaries would be more efficient than the current theory 
predicts. Our revision of the torque density prescription 
should have important effect on the self-consistent cal- 
culation of the gap density profile and, consequently, on 
the orbital evolution of the perturber (Type II migration 
speed for planets and SMBH inspiral rate). Future ex- 
tension of our results to the full cylindrical geometry will 
allow direct applications of the non-uniform disk theory 
to systems with wide gaps or cavities. 
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APPENDIX 
SHIFTED LINDLBLAD RESONANCES IN THE SHEARING-SHEET APPROXIMATION 
The angular frequency of the non-uniform isothermal disk in full cylindrical geometry is given by 



Vi\v) = n 2 K (r) + 



<91n£n 

dr 



(Al) 
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Fig. 9. — Torque density outside the gap for three different models J23D with S m i n = 3 X 1 0~ 4 , and different n and S as specified in each 
panel. Black curves show our numerical results while the red dashed curves show the fit l)BH l with the input parameter a indicated in each 
panel. 



where Q,k in the regular Keplerian frequency. The modified epicyclic frequency is then given by 



k 2 = 4fi 2 + 2rO— - = n 2 K (r) + c 



,<9 2 ln£ 3c 2 <91nE 



dr s dr 2 r dr 

2 



(A2) 



The locations r^ at which the Lindblad resonance condition m 2 [0(ri,) — fl p ] = « 2 (ri,) (here fi p = £Ik{t p ) is the 
pattern speed of the gravitational perturbation and m = k y r p is the wavenumber) are given by the following relation: 



^ 2 K(r L ) + 



c 2 d In S 



di 



i Lr 



2 f n , 2/'5 2 lnE , 3 51nE 



dr 2 



Oi 



(A3) 



Expressing m via k y we find the following relation between k y and corresponding 77,, which does not make any 
approximations and is valid in full cylindrical geometry: 



k y {r L ) = ± 



«(r L ) 

i IpTp 




h 2 dlnZr 



r dr 



- 1 



, K 2 {r L ) = n 



r„y | ft2 / 9 2 lnS | 3 91nS 

7'^ y V 9r 2 r 9r 



(A4) 



In the shearing sheet approximation one takes rj, = r p + Xl, expands r p jrj 1 to linear order in xl/t p and then takes 
the limit r p — > 00. As a result one recovers the resonance condition in the form (|24|) with k given by equation ([25]>. 

A ONE-PARAMETER TORQUE DENSITY FIT 

Given the numerical challenges involved in the calculation of the torque density in a general non-uniform disk we 
try to provide an empirical fit for dT/dx. Our procedure is based on the following logic. 

For every gap model considered in this work we find that the torque density per unit surface density £0 decays 
exponentially outside the gap until it becomes negative because of the negative torque density phenomenon (RP12). 
In particular, for our fiducial model the fall-off is described reasonably well by exp(— 0.8x/h) for x/h > 5 as shown 
in Figure [6l Also, as discussed in §6. 2\ [Tithe full integrated torque T(oo) coincides within ~ 10% with the prediction 
based on simple GT80 prescription (0| and (|3"2"j) . We use this property to get a normalization to the exponential fall-off 
outside the gap. Putting these two ideas together we come up with the following fit for the torque density outside the 
gap (x > 5): 



— = C(a)sign(x 
dx 






n 2 h 4 



f™ Z Q (x)e- ax / h dx ' 



(Bl) 



where a is a dimensionless free parameter that controls the exponential fall-off and takes value in the range ~ 0.5 — 1 
for all our numerical calculations as we show next. 

In Figure @, we show how this fit performs for different gap profiles as we vary the width or steepness of the density 
profiles (the torque density outside the gap does not depend on its density contrast S m i n , see £|7.3[) . In the left panel, 
we consider our fiducial model and the fitting formula using a = 0.8. Our prescription fits the numerical result pretty 
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well from bh < x < 15h, while beyond 15ft. the torque density becomes negative and the correct asymptotic behavior 
is given by the analytical expression — 0.63T, (x)/x 4 , see RP12 (in proper units). 

In the middle panel, we increase the width of the gap by taking 5 = 7h and show a fit with a — 0.58, which works 
well for 7h < x < 22h, before torque density changes sign. 

In the right panel, we increase the steepness of the gap profile by using n — 4 and show the best fit (|B1[) . which now 
requires a = 0.95. In this case, the fit is worse than in the previous examples with relative errors up to < 20% at the 
peak of the torque density, but it still does a better job than the asymptotic expression from GT 80. 

Unfortunately, at the moment we cannot predict a value of the parameter a featured in equation (|B1[) for a particular 
gap profile. But the results shown in Figure [9] suggest that the torque density decay is faster and the fits demands 
higher values of a for boxier density profiles (i.e. for higher n) and for narrower gaps (i.e. for lower S). 



